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Abstract 

Novel electronic states resulting from entangled spin and orbital degrees of freedom are hallmarks 
of strongly correlated /-electron systems. A spectacular example is the so-called 'hidden-order' 
phase transition in the heavy-electron metal URu2Si2, which is characterized by the huge amount 
of entropy lost at Tho = 17.5 K 12 . However, no evidence of magnetic/structural phase transi- 
tion has been found below T HO so far. The origin of the hidden-order phase transition has been 
a long-standing mystery in condensed matter physics. Here, based on a first-principles theoretical 
approach, we examine the complete set of multipole correlations allowed in this material. The re- 
sults uncover that the hidden-order parameter is a rank-5 multipole (dotriacontapole) order with 
'nematic' E~ symmetry, which exhibits staggered pseudospin moments along the [110] direction. 
This naturally provides comprehensive explanations of all key features in the hidden-order phase 
including anisotropic magnetic excitations, nearly degenerate antiferromagnetic-ordered state, and 
spontaneous rotational-symmetry breaking. 
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In the rare-earth and actinide compounds, /-electrons behave like well-localized moments at 
high temperatures. As the temperature is lowered, /-electrons begin to delocalize due to the 
hybridization with conduction electron wavefunctions. At yet lower temperatures the /-electrons 
become itinerant, forming a narrow conduction band with heavy effective electron mass, which 
is largely enhanced from the free-electron mass. Notable many-body effects within the narrow 
band lead to a plethora of fascinating physical phenomena including multipole order, quantum 
phase transition and unconventional superconductivity. Among them, perhaps the appearance of 
a 'hidden-order' (HO) state in URu 2 Si 2 is one of the most mysterious phenomena. Identification 
of the microscopic order parameter and mechanism that derives the HO transition continue to be a 
central question in the strongly correlated /-electron systems. 

There are several unique features that appear to be clues for understanding the HO in 
URu 2 Si 2 . In the paramagnetic state above T HO , the magnetic susceptibility exhibits the Ising-like 
anisotropy 1 ' 3 . In the HO state below T H o, an electronic excitation gap is formed on a large portion 
of the Fermi surface (FS) 4 5 and most of the carriers disappear 6-8 . Closely related to this, the gap 
formation also occurs in the magnetic excitation spectra at commensurate and incommensurate 
wave numbers, Qc = (0 1) and Qic = (0.6 0), respectively, as revealed by the neutron inelastic 
scattering 911 . The HO ground state changes to the large-moment antiferromagnetic (AFM) state 
with the ordering vector Q c upon applying hydrostatic pressure 1214 , but the FS has a striking 
similarity between these different phases 416 , implying that the HO is nearly degenerate with the 
AFM order. The magnetic torque measurements reveal the 'nematicity', which breaks the in-plane 
rotational (tetragonal) symmetry in the HO 17 . The challenge for the theory has been to identify the 
order parameter which explains all the above key features. 

The theories that have been proposed to describe the HO state can be divided into two pre- 
vailing approaches; one is based on the localized 5/-electron model and the other the itinerant 
one 18 . Recent angle-resolved-photoemission-spectroscopy results clearly demonstrate that all 5/ 
electrons are itinerant 19 and the crystalline electric field, which is a signature of the localized na- 
ture, has never been observed. Moreover, the nuclear magnetic resonance (NMR) measurements 20 
show a formation of the coherent heavy-electron state well above T HO - Therefore it is natural 
to discuss the electronic structure based on the itinerant picture. However, reliable calculation of 
the physical quantities by taking into account the complicated band structure is a difficult task. 
For this purpose, we use a state-of-the-art ab initio downfolding 21 and dissect the electronic struc- 
ture obtained from the density-functional theory (DFT) calculations. The obtained tight-binding 

2 



Hamiltonian is constructed from 56 orbitals of U 5/, U Qd, Ru Ad and Si 3p. Introducing the 
on-site Coulomb interactions between 5/-electrons, we obtain a realistic itinerant model, i.e. 56- 
band Anderson-lattice model including the spin-orbit interaction. Based on this realistic model 
Hamiltonian, magnetic and multipole correlations are analyzed by the random-phase approxima- 
tion (RPA) and beyond. To account for the mass renormalization effect in the Fermi-liquid theory, 
the energy and temperature scale is reduced by a factor of 10 1,4,16 throughout this study, which 
makes comparisons to the experiments straightforward. 

Figure 1 displays the paramagnetic FS and the band structure near the Fermi level, respectively. 
The energy bands crossing the Fermi level have mainly the total angular momentum j = 5/2 
multiplet of U 5/. Each j z component of j = 5/2 multiplet is weighed by color. It turns out that 
each separated FS is mainly composed of a rather specific j z component without large mixing, 
except for the outer FS around Z point (Fig. 1). Such a j z component map is quite useful in that 
we are able to capture valuable information such as which parts of the FS play an essential role for 
the HO formation. Indeed, the disentanglement of FS orbital characters has also been an important 
theoretical advance to understand the electronic properties in iron-pnictide superconductors 22 . 

First we discuss the RPA analysis of rank-1 (dipole) correlation, which is the conventional 
static magnetic correlation. The regime with j z = ±5/2 shown by red in outer FS around Z 
point is well nested with outer FS around T point by the vector Q c as indicated by the arrow in 
Fig. I 3 . This nesting gives rise to a sharp peak of the correlation parallel to the c axis (dipole 
J z ) at Z (0 1) shown in Fig. 2a. Another salient feature is the hump structure at around (0.6 
0) and the equivalent points, whose Q-vectors coincide with Qic- We point out that these peak 
and hump structures in the paramagnetic phase are directly related to the magnetic excitation gap 
at Q c and Q IC in the HO phase 11 because the gap opening occurs at the nested parts of the FS. 
We also note that the in-plane magnetic correlations, J x and J y , are not enhanced in contrast to 
J z (Fig. 2b), which is also consistent with the Ising-like magnetic susceptibility 1,3 and polarized- 
neutron measurements 10 . 

Next we examine the higher-rank multipole correlations. According to the group theory, there 
are 36 multipole moments up to the fifth rank in j = 5/2 subspace (Table SI). Figures 2 c-f depict 
the correlations between the basis functions belonging to rank-2 (quadrupole), 3 (octupole), 4 (hex- 
adecapole) and 5 (dotriacontapole), respectively. What is remarkable is that similar to the dipole 
J z , the Qc correlation at Z-point is strongly enhanced in some cases such as O 2 o (rank-2), H x ^ b 
(rank-4), D 4 (rank-5), etc. Generally, these bases are mixed in the tetragonal symmetry, as shown 
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by finite off-diagonal terms (red lines). The multipole correlations obtained by the diagonalization 
is depicted in Fig. 2g, in which each correlation at Q c is classified by the irreducible represen- 
tations and the dominant component is denoted in parentheses. At low temperatures, A 2 (J Z ), 
E~(D x { y \) and A{{D A ) symmetries exhibit the first, second and third strongest enhancement. The 
former and the latter two correspond to the AFM and dotriacontapole states, respectively. Within 
the RPA, the AFM state always overcomes the dotriacontapole states. To go beyond the RPA, we 
take into account the mode-mode coupling by including the Maki-Thompson type vertex correc- 
tions. For this purpose, we calculate the maximum eigenvalue A of the Bethe-Salpeter equation 
for a staggered particle-hole pairing with use of the RPA result as the kernel. At A = 1, a phase 
transition occurs from the paramagnetic to the corresponding eigenstate. Temperature dependence 
of A of each symmetry is shown in Fig. 3a. As the temperature is lowered, A of E~ (D x < y \) is most 
strongly enhanced and the condition A = 1 is fulfilled at finite temperature, indicating a phase 
transition to the E~{D X ^) state. We emphasize that E~{D X ^) symmetry breaks the in-plane 
fourfold symmetry, which naturally accounts for the 'nematicity' observed in the magnetic torque 
results 17 ' 24 . In addition, E~(D x ( y )) state breaks the time reversal symmetry, which is consistent 
with the NMR measurements 20 . These lead us to conclude that the HO is E~(D X ^) dotriacon- 
tapole order. 

The present calculations also reproduce well other key features of the HO, i.e. near degeneracy 
of the HO and AFM states and the anisotropic temperature dependence of the uniform suscepti- 
bility. Fig. 3a demonstrates that A of E~{D X ^) is very close to that of A 2 (J Z ). This indicates 
that both states are nearly degenerate and small perturbation can change the HO to AFM state. 
Indeed, we can construct a phase diagram by tuning the interactions 21 , which is consistent with 
the pressure-temperature phase diagram (Fig. 3b). The temperature dependence of the uniform 
susceptibility x c (0) parallel to the c axis exhibits a broad maximum at around 40K, whilst Xab(0) 
perpendicular to the c axis is smaller and nearly temperature independent (Fig. 3c), in good agree- 
ment with experiments 13 . The low-temperature decrease of x c (0) arises from the deep dip struc- 
ture in the density of states (DOS) near the Fermi level (Fig. Sib). The Ising-like susceptibility 
including its temperature dependence has been discussed in terms of the crystalline electric field 
excitations of the localized 5 /-electrons so far. However, the present results demonstrate that the 
susceptibility can be well accounted for by the itinerant scenario. 

Why is such a high rank multipole state (rank-5) with E~(D x ( y ))-symmetry realized in 
URu 2 Si 2 ? In the paramagnetic state, the FS nesting with Q c vector plays an essential role on 
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the multipole fluctuations. What is crucially important is that the FS regions connected by this 
Qc are dominated by the ±5/2 components, as shown in Fig. 1. In this situation, we can con- 
sider a subspace consisting of only two components j z = 5/2 and —5/2, which allows us to map 
j z = ±5/2 to pseudospin | and |. Then the dipole J z is described by the Pauli matrix a z , as it has 
only diagonal elements corresponding to ±5/2 <^> ±5/2. In contrast, D x ^ is given by cr x ( y ), rep- 
resenting off-diagonal components describing the ±5/2 <^> =f5/2 transition 21 , which accompanies 
the angular momentum change of 5h allowed only in rank-5. In this pseudospin space, the stag- 
gered J z state corresponds to the Neel order along the c axis. On the other hand, the D x ^ state 
corresponds to the in-plane order breaking the rotational symmetry, where in-plane pseudospin 
moments are antiferromagnetically coupled along the c axis (Fig. 3d). Thus the pressure-induced 
first-order transition from HO to AFM state can be explained by the pseudospin staggered mo- 
ment flip from the in-plane to the out-of-plane direction. The experimentally observed 'nematic- 
ity' along the [1 10] -direction corresponds to the linear combination -D[n ] = ^75 (Ac + D y ) of the 
two-fold degenerate D x and D y . The HO parameter is then represented by 

0[iio](fc)= (fL a [uo]h+Qcfi), °lno] = — + J? V , (1) 

where fk a is an annihilation operator for an /-electron with momentum k and pseudospin a. It 
should be noted that under in-plane 180° rotation, the pseudospins change their direction, which 
discriminates this state from a nematic phase in the strict sense. However, its staggered nature 
leads to the twofold 'nematic' symmetry of the bulk susceptibility as observed experimentally. 

Figure 4a displays the FS in the HO state, which is calculated by applying the effective multi- 
pole field so as to open the gap of 4 meV observed by the scanning tunneling microscopy 4 ' 5 . The 
lattice doubling in the AFM phase with Q c also occurs in the HO phase. Most part of the FS 
having j z = ±5/2 components disappears as a result of the gap opening at the nested parts of 
the paramagnetic FS. Around the T-point, a small electron and a large hole (a) pockets, the FS 
with a cage-like structure and four electron pockets (j3) exist. The FS in the HO phase bears a 
striking resemblance to that in the AFM state, consistent with the quantum oscillation measure- 
ments. However, the broken fourfold symmetry in the HO state can be seen clearly in the FS with 
cage-like structure (Fig. 4b), in sharp contrast to the AFM state 21 . 

The present approach based on the first-principles calculation is able to give a comprehensive 
explanation to the problem of HO, which has been a quarter century mystery. Why has the HO 
been hidden for a long time? The reason appears to be that the order parameter of the present 
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high-rank multipole state is extremely difficult to detect directly by the conventional experimental 
techniques, such as resonant X-ray and neutron measurements. The itinerant multipole ordering 
with 'nematicity' revealed in the present study is a new type of electron ordering, which is expected 
to be ubiquitously present in the strongly correlated electron systems 25 when spin and orbital 
degrees of freedom are entangled. 
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FIG. 1. Paramagnetic FS and energy band dispersion colored by weight of j z component. Red, green 
and blue color gauges correspond to j z = ±5/2, ±3/2 and ±1/2 components, respectively. The FS is 
constructed from two hole FSs around Z, and the other four electron FSs. Small (blue) electron pockets 
centered at X and T are constructed from j z = ±1/2 component, and the inner (green) hole pocket around 
Z is from j z = ±3/2. The outer hole pocket around Z is a hybridized band between ±3/2 and ±5/2. The 
outer electron FS around V is mainly composed of j z = ±5/2, and partially hybridized with ±1/2. Two 
outer FSs around V and Z are partially nested with Qc = (0 1) indicated by arrow. 

FIG. 2. (a) Momentum dependence of magnetic susceptibility at T = 12K for J z . A complete set of 
multipole correlations at T = 12K are shown along high-symmetry line for rank-1 dipole (b), rank-2 
quadrupole (c), rank-3 octupole (d), rank-4 hexadecapole (e) and rank-5 dotriacontapole (f) basis functions. 
Off-diagonal correlations between different bases are also shown in red curves, (g) Diagonalized multipole 
correlations as a function of temperatures, where D x r y \ = (D x / y \ al + D x i y \ a2 + D x / y \ b )/y/3. These all 
correlations have been obtained within the RPA calculations for U = U' ~ 2.3 and J = J' = in units of 
l/pf, where pf is the total /-electron DOS at the Fermi level. 

FIG. 3. (a) Temperature dependence of the maximum eigenvalue A of the Bethe-Salpeter equation beyond 
the RPA including the Maki-Thompson type vertex corrections for U = U' ~ 2.4 and J = J' = 0. Note 
that A = 1 gives a phase transition temperature to the corresponding eigenstate. (b) By increasing the 
Hund's coupling J 21 , the small difference of A between the rank-5 E~(D x t y \) and rank-1 A 2 (J z ) states 
can be reversed, which may account for the pressure-induced AFM state, (c) Temperature dependence of 
uniform susceptibilities parallel and perpendicular to the c axis calculated beyond the RPA. (d) Schematic 
configurations of the ±5/2 pseudospin moments in the HO (left) and AFM (right) states are shown by the 
arrows. In both states, the pseudospins order antiferromagnetically along the c axis, but the direction of the 
staggered moments in the pseudospin space differs between the two: along [1 10] for the HO state and along 
[001] for the AFM state (center). 

FIG. 4. £ ,_ (Dr 110 i)-state FS (a) and the band dispersion (b) colored by weight of the j z components, (a) 

The Bulilluoin zone is folded with Qc = (0 1) and a two-dimensional cut of FSs is shown in the k z = 

plane. The electronic structure exhibits the in-plane four-fold symmetry breaking, which can be most easily 

seen in the cage FS. Two FSs around M (blue and green lines in the right bottom) have almost no splitting 

along XM, which is in contrast to the large splitting found for the AT, AFM and the E + state 21 . (b) The 

dispersion along TM line (left panel) highlights the excitation gap of ~ 4meV (arrow). The enlarged figure 
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near the cage FS (right panel) shows a pronounced anisotropic TM dispersion between E (red) and £' 



SUPPLEMENTARY INFORMATION 



Emergent Rank-5 'Nematic' Order in URu 2 Si 2 



I. ELECTRONIC BAND-STRUCTURE CALCULATION AND AB INITIO DOWNFOLDING 

First, we perform the ab initio band- structure calculation in the paramagnetic state of URu 2 Si 2 
using the wien2k package 1 , in which the relativistic full-potential (linearized) augmented plane- 
wave (FLAPW) + local orbitals method is implemented. The crystallographical parameters are 
the space group No. 139, IA/mmm, the lattice constants, a = 4.126A, c = 9.568A, and Si internal 
position, z = 0.371 2 . Red line in Fig. Sla depicts the electronic band structure, which is very 
similar to previous works 3 4 . The band structure is composed of U 5/, U 6d, Ru 4c? and Si 3p 
orbitals. U 5/ band is located where the bottom of U Qd band overlap with the top of Ru Ad 
band. On the basis of the LS basis set, spanned by these orbitals and spin (| and I), we carry out 
WANNIER90 code 5 via wien2wannier interface 6 . With one-shot calculation of the wannier90, 
we obtain a real-space representation of Kohn-Sham equation, i.e., the tight-binding Hamiltonian 
in 14 + 10 + 10 x 2 + 6 x 2 = 56 orbital bases. We do not make maximally-localized Wannier 
functions (MLWFs) to preserve the on-site symmetry of the p, d, and / orbitals, since construction 
of MLWFs leads to a mixing between up and down spin components at each MLWF The obtained 
Wannier fit (blue dashed line) is well consistent with the original band structure. The tight-binding 
Hamiltonian is written as 



where f^e (cki) is an annihilation operator for an /- (conduction-) electron with momentum k 
and orbital i. Superscripts / and cond. in the sums denote all 14 /-orbitals and 42 conduction 
bands, respectively. The average Wannier spread of each orbital is 0.9, 2.9, 1.6, 3.0 in the atomic 
unit for U 5/, U 6d, Ru Ad, Si 3p. As expected, U 5/-orbital wavefunctions are well confined 
at an Uranium site. This implies that the on-site Coulomb interactions between /-electrons is the 
largest and the most important. The /-orbitals in the LS basis are transformed into the J basis with 
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FIG. SI. Band structure (a) along high-symmetry line. Red line is the result of DFT calculation by WIEN2K. 
Blue dashed line is the Wannier fit. The dispersion below ~ 2eV is reproduced completely. The J-resolved 
DOS (b). The left (right) down figure depicts the partial DOS in the j = 5/2 (7/2) manifold. The Fermi 
level is located in a dip structure in j = 5/2 bands, j = 7/2 bands lies leV higher. The Fermi surface 
colored by the Fermi velocity (c). The Fermi velocity is large around k z axis. 
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where subscripts £ and j denote orbitals in the LS and J bases, respectively. The LS basis set is 
given by the direct product of the orbital space, {z 3 , xz 2 ,yz 2 , z(x 2 — y 2 ), xyz, x(x 2 — 3y 2 ), y(3x 2 — 
y 2 )}, and the spin space (|, 1). The J basis set is the direct sum of the total angular momentum 
j = 5/2 space, (+5/2,+3/2,+l/2,-l/2,-3/2,-5/2), and the j = 7/2 space, (+7/2,+5/2,+3/2,+l/2,- 
1/2,-3/2,-5/2,-7/2). 

Figures Sib shows the J-resolved density of states (DOS). The DOS near the Fermi level, 
Ep = OeV, is dominated by the j = 5/2 multiplet. The weight of j = 7/2 orbitals exists around 
leV higher due to the local crystalline electric field and the strong spin-orbit interaction, XL ■ S 
(A = 0.24eV). It should be noted that the Fermi level lies at a deep dip in the DOS. As shown latter, 
this feature is important for the temperature dependence of the uniform magnetic susceptibility. 
The occupation numbers for these j = 5/2 and 7/2 multiplets are 2.07 and 0.64, respectively. The 
orbital occupancy in the j = 5/2 multiplet is 0.72, 0.86, and 0.49 for J z = ±5/2, ±3/2, and 
±1/2, respectively. 

Figure Sic is the Fermi surface (FS) colored by the Fermi velocity. Blue color means heavy 
band mass. The FS around k z axis possesses light band mass due to a large mixing with Ru Ad 
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bands. 

Finally, we construct a realistic Anderson lattice model by adding the on-site Coulomb in- 
teractions between /-electrons to this ab initio tight-binding Hamiltonian, H . The interaction 
parameters (U, U', J, J') are introduced in a conventional form in the LS basis. 

H ' = ^J2Y1 fitafmfi^fi^ ( S2a ) 

it a 

+ y J] ^ fuafima'fima'fita (S2b) 
ity^m aa' 

+ ^ X^ X^ fie a fj ma i flier' fima (S2c) 
iiy^m era' 

J' 

+ y fi£ a fi£afimafima, (S2d) 

itj^m a 

where er = ± and a = — er. We measure these interaction parameters in units of 1/p/, where 
P/ = 4.20 (States/eV) is the total f-electron DOS at E F . Hereafter, to account for the mass 
renormalization effect, the energy and temperature scale is reduced by a factor of 10. 



H. GROUP-THEORETICAL ARGUMENT 

We here define one-particle operators for multipole moments from the Group-theoretical ar- 
gument. Since the DOS near the Fermi level is dominated by the j = 5/2 multiplet, we neglect 
the higher- level j = 7/2 multiplet, and consider only the j = 5/2 manifold. In the six states 
(j z = ±5/2, ±3/2, ±1/2), the irreducible tensors are completely available up to rank 5; dipole 
(rank 1), quadrupole (rank 2), octupole (rank 3), hexadecapole (rank 4), and dotriacontapole (rank 
5) moments. 

(k) 

The irreducible tensor for rank k has Ik + 1 components Jg which meet the following rela- 
tions, 

Jz,J { q k) ] =qJ q [k \ (S3) 

J ± , jf ] = y/(kTq)(k±q+l)J% (S4) 

where J± = J x ±iJ y are raising and lowering operators, and J M is /j, component of the total angular 
momentum operator. In Table S 1 , we show multipole moments and Hermite bases belonging to 
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TABLE SI. Definition of multipole moments under the tetragonal D^h symmetry in the j = 5/2 subspace, 
which are completely described with 36 bases up to rank 5; dipole (J, rank 1), quadrupole (O, rank 2), 
octupole (T, rank 3), hexadecapole (H, rank 4), and dotriacontapole (D, rank 5), including monopole 
(rank 0). These bases are normalized by Eq.(S21). Superscripts ± of irreducible representations denote the 
parity under the time-reversal transformation. D x ^ = (D x ^ al + D x ^ a2 + D x ^ b ) / \/3. 



Symmetry Notation 


Basis 




o 20 


J o 




H 






H 4 


(V5J^-V7J^ + )/VT2 


A% 


H za 


J [4]- 


Bt 


O22 


f(2) 
J [2} + 




H 2 


J [2} + 


Bt 


Qxy 


f(2) 
J [2\- 




H z b 


j(4) 
J [2]- 


E+ 


Oyz ) O zx 


f(2) f(2) 




H X aj By a 


/ z(A\ /zr z( A\ \ 1 1 — / z(A\ 1 — z( A\ \ 1 1 — 

(41+ + ^ J S+)/^' -(41- - ^41-)/^ 




H x b, Hyb 




A', 


D A 


J [i]- 


A- 2 


Jz 


7(1) 

J 




T Z a 


7(3) 
J 




Dzal 


7(5) 
J o 




D Z a2 


7(0) 
J [4] + 


BT 


q^xyz 

D 2 


f(3) 

J m- 
7(5) 
J m- 


B2 




7(3) 

J [2] + 




D zb 


f(5) 
J [2] + 


E~ 


Jxi Jy 


J [l]-> + 




T X a 1 Tya 


(^4f_ - Vsjf^J/A -(v/sjgV + Vajg^/Vs 




Txbi Tyb 


-(^4f_ + VsjgU/Vs, -(^4f + - Vsjg^/Vs 




D X al j D ya i 


(3^J ( ( |_ - V70J§_ + 2 % /T5jg ) _)/16, (3VT44V + V70J§+ + 2 v /l5jf 1 , j ) + )/16 




D X a2 , D ya 2 


(VWJ§_ + 9y/2 + 2^4?-)/ 16 ' (^41+ " 9 ^4l+ + 2 ^ IJ [i]+)/ 16 




D x bi Dyb 


(V30Jg_ + y/6J$_ - 2y/7J®_ )/8, {V30jg + - y/6J§} + - 2y/7J$ + )/8 
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each irreducible representation of the tetragonal symmetry (D ih )' '. Here J^ ± is defined by 



J [2p^l] + 



J, 



(fc) 

[2p-l]- 
?(*) _ 

t(*) . 



1 f_ t(*) , r(fc) 
* f _ _|_ 7W ^ 



(S5a) 
(S5b) 
(S5c) 
(S5d) 



for positive integer, p. Applying the Wigner-Eckart theorem, we obtain the matrix elements of 



t(*0 

Jq ' as 



(S6) 



where (jWJ^Wj) is the reduced matrix element, and (jj z \jj' z kq) is the Clebsch-Gordan coeffi- 
cients. These are easily evaluated in terms of Mathematica. Thus we obtain numerically all 
representation matrices in Table SI. For example, matrix elements of higher-rank multipoles, H xb , 
D 4 , and D x(y) = (D x{y)al + D x(y)a2 + D x(y)b )/^/3, are explicitly given by 



H 



xb 



( 0.09i -0.47i O.OoA 

-0.09i — 0.15z O.OOi 

0.15i O.OOi 0.47i 

0.47i O.OOi 0.15i 

O.OOi -0.15i -0.09i 

\ O.OOi -0.47i 0.09i / 



(S7) 



Da 



( -0.50i N 

O.OOi 0.50i 
O.OOi O.OOi 
O.OOi O.OOi 
0.50i O.OOi 
\ -0.50i / 



(S8) 
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D x = 



( 0.02 0.11 0.65\ 

0.02 -0.08 -0.08 

-0.08 0.11 0.11 

0.11 0.11 -0.08 

-0.18 -0.08 0.02 

\0.65 0.11 0.02 / 

/ -0.02i O.lli -0.65i\ 
0.02i 0.08i -0.08i 

-0.08i — O.lli O.lli 

-O.lli O.lli 0.08i 

0.18i -0.08i -0.02i 

\ 0.65i -O.lli 0.02i J 
where those norms are normalized by Eq.(S21). 



(S9) 



(S10) 



III. RPA AND MAGNETIC CORRELATIONS 



A. Formalism 



First, we calculate one-particle Green functions in the LS basis, 



G em (k,iu n ) = -((fktf; 



km 



= -["dre^(T r [f k e(r)ft m (0)]) 
Jo 



(SI la) 
(SI lb) 



where £ and m denote both /-orbital and spin quantum number. The non-interacting static full 
susceptibility is given by 



k,n 



lU n ). 



(S12) 



This is written in the 14 2 x 14 2 matrix form, x°(q), with the {£m} row and the {£'m'} column. In 
this case, the RPA susceptibility is given by 



i-x°(?)r° " X°(q), 



(SI 3a) 
(S13b) 
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where T is the bare interactions between /-electrons, 



££,mm 


~~ 1 £+7£+7,m+7m+7 ~ 


J-U', 


r° 

£m,£m 


— 1 £+7m+7,e+7m+7 ~ 


U' - J, 


L U,m+7m+7 


1 £+7£+7,mm 


-u>, 


1 £m,m+7£+7 


- r° 

- 1 - £+7m+7,m£ 


-J', 


1 £m,£+7m+7 


- r° 

~~ - 1 - £+7ra+7,fon 


-J, 


1 ££+7,mm+7 


- r° 

x £+7£,m+7m 


J, 


r° 

L £m+7,m£+7 


- r° 

- 1 - £+7m,m+7£ 


J', 


r° 

L £m+7,£m+7 


_ r o 

- 1 - £+7ra,£+7m 


U', 


1 £,£,£+7,£+7 


— 1 £+7,£+7,£,£ 


-u, 


1 £,£+7,£,£+7 


— 1 £+7,£,£+7,£ 


u, 



with orbital indexes £ ^ m G 1 ~ 7. Correlation between multipole A and multipole B is 
evaluated by the product of the above susceptibility and the representation matrices, 



{(A,B))= A me xZ A e m/ (q)B e/m , (S14) 

££'mm' 

In the LS basis set, the spin moment with /i = x, y, z is the product of the unit matrix in the 
orbital space and the Pauli matrices in the spin space, 1 ® The orbital moment is the product 
of the unit matrix in the spin space and L M ; 



/ VQi 


— y/6i 



\ 



o -^pi 



2 1 







2 4 



V6a 



(S15) 
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(S16) 
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L z = 
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(S17) 
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-2i 
2i 

-Si 
3i 

The total magnetic moment M M is defined by L M + 25^. 

Next, let us consider only the j = 5/2 subspace, neglecting the j = 7/2 subspace. The 
magnetic moment is approximated as M M = gj^, where g = 6/7 is the Lande g-factor. The 
magnetic dipole moment, J M , in the j = 5/2 subspace is given by 



u 2 

^ 

y/2 



\ 



3 
2 



V 



V2 


f o J 



x/2 ^ 



(SI 8) 
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2 * 



\/2i 
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-§i 




§i 


-\/2i 
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(S20) 
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In actual calculations, since we have already obtained the RPA susceptibility in the LS bases, we 
need to transfer these J-base representation matrices into the LS-basis ones, J2jj> c tj^/ c mj' by 
the unitary matrix cej. 



B. Magnetic correlations 

We carry out the RPA analysis for several parameters in 32 x 32 x 8 /c-meshes. For simplicity, 
we show in the present study only some results for U — U' ~ 2.3 and J = J' = 0. The obtained 
characteristic features are barely changed for J = J' > 0. Mainly, the RPA treatment remarkably 
enhances the magnetic character in the non-interacting system. Figure S2 is the magnetic cor- 
relations along high-symmetry line. Figure S2a demonstrates a remarkable magnetic anisotropy, 
namely, the L z correlation ((L Z ,L Z )) larger than the L x ^ correlations ((L^), L^)}). The L z 
correlation indicates a peak structure at Z (1 0) point and a hump structure at (0.6 0), while 
the L x ty\ correlation is featureless. In these correlations, the correlation between orbital moments 
(red line) are much larger than spin-spin correlations (green line). The correlations between spin 
and orbital moments (blue line) are negative, and compensate the large orbital-orbital correlations. 
This means that the large orbital moment is compensated by the anti-parallel spin moment. Gen- 
erally, this is the case in the Uranium compounds, since the electron occupation of / orbitals is 
less than half. In Fig. S2b, we compare the total magnetic correlations with the gj correlations 
confined in the j = 5/2 subspace. As expected from the fact that the DOS near E F is dominated 
by the j = 5/2 components, the total magnetic correlations are overall explained by the j = 5/2 
correlations. Thus we expect that we can neglect the effect of the j = 7/2 subspace also in the 
higher-rank multipoles. 

Figure S3 depicts temperature dependence of the uniform magnetic correlations within the 
RPA. The J z correlation x^ PA (0) parallel to the c axis is larger than the in-plane J x ^ correlation 
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FIG. S2. LS-resolved magnetic correlations (a) along high-symmetry line, where M = L + 25. The orbital 
contribution is much larger than spin contribution (green line). Correlation between orbital and spin (blue 
line) is negative due to the less than half, j = 5/2 subspace contribution (b) to the total magnetic correla- 
tions along high-symmetry line. Differences between ((M^, M^)) correlation and g 2 (^J^, J^)) correlation 
come from the contribution of j = 7/2 manifold, which provides only small constant shift. 




1 1 1 1 1 1 

100 200 300 

T(K) 

FIG. S3. Temperature dependence of the uniform susceptibility within the RPA. The J z correlation parallel 
to the c axis (red), x^ PA (0), is larger than the in-plane J x t y ) correlation (blue), Xa PA (0)- B° tri correlations 
show a broad maximum at around 50K. 

Xafe PA (0)- Both correlations show a broad maximum at around 50K. Such temperature dependence 
is similar to the experimental result for the c-axis magnetic susceptibility. This feature arises from 
the fact that the DOS possesses a deep dip structure near E F . In calculations beyond the RPA 
(Sec. V), the uniform susceptibility Xc(0) parallel to the c axis is enhanced, while the in-plane 
Xab(0) is suppressed. In this case we obtain more remarkable Ising-like anisotropy including the 
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temperature-independent behavior in Xab(O) (Fig. 3c). 



IV. MULTIPOLE CORRELATIONS 

We obtain higher-rank multipole correlations (Fig. 2) in the same way as the J M correlations. 
In order to compare different-rank multipole correlations, we here normalize those norms as 



where Q fm represents matrix elements of multipole moments in Table SI. Owing to this normal- 
ization, magnitude of dipole correlations in Fig. 2a becomes 2/(35g 2 ) smaller than in Fig. S2b. 
Figure 2 shows the following features. In quadrupole correlations, O 2 o and O yz ( zx ) are dominant 
at Z point. Octupole correlations have featureless Q dependence. In hexadecapole correlations, 
H x { y )p is dominant. In dotriacontapole correlations, D±, D za2 , D x r y \ alj2 an d D x i y \p show a peak 
structure at Z point. It should be noted that the off-diagonal part (red lines) is large in dotria- 
contapole correlations. This means that different bases have a large mixing. Such mixing is not 
restricted in the same rank. Thus we need to diagonalize the full susceptibility Xtm,t'm'{Qc) itself 
in the 14 2 x 14 2 matrix form. The results were shown in Fig. 2g. The obtained nth eigenvector, 
\I/" m , can be expanded in terms of 36 multipoles, Q\ m , shown in Table SI; 




(S21) 




(S22) 



Each component of A™ was obtained as the following, 




0.109 Jj, - 0MlT ya + 0.190Tj, 6 
+0M2D yal + 0A13D ya2 + 0.540 A, 6 , 



0.824J Z + 0.281T za - 0.298D zal + 0.391D. 




0.109,4 - 0.3417; a + 0.190T, 



+0.612D xa i + 0.413^2 + 0.540LU 



(S23) 




-0.266O 



L> 4 , 



yz 



0.372H xa + 0.889if x6 



4' 
E+ 



-0.2QQO ZX - 0.372H ya + 0.889H yb , 
-0.2910^ + 0.957H zb . 
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2(2') 4(4'') 

E~ and E + representations are twofold degenerate, \[> and . The dominant order pa- 
rameter is dipole, E~ dotriacontapole, A[ dotriacontapole, and E + hexadecapole, and 
hexadecapole, in order. The E~ state includes comparably D x i y \ a i, D x i y \ a2 , and D x i y \ b . Thus we 
have introduced a new bases D x ^ = (D x ^ al + D x ^ a2 + Ac(y)&) / in the present study. 

V. BEYOND RPA AND PSUDOSPIN REPRESENTATIONS 
A. beyond RPA 

Generally, RPA enhances a magnetic channel, but depresses a charge channel. This trend is 
improved by including higher-order fluctuations beyond the RPA, i.e., the mode-mode coupling 
terms such as Maki-Thompson type vertex corrections. In multi-orbital systems, such mode-mode 
coupling term mix also many different channels. This implies that a channel which is not en- 
hanced within the RPA may develop. Therefore calculations beyond the RPA is very important 
especially in the multi-orbital systems. We here examine Maki-Thompson type mode-mode cou- 
pling effect to search a possibility that some high-rank multipole correlation overcomes the dipole 
J z correlation, which is dominantly enhanced within the RPA. 

The present Maki-Thompson type diagrams are shown in Fig. S4a. The wavy line represents 
an effective interaction including the RPA results, 

Vem,e'm'(<l) — r°m,*'m' + ^m,£ 1 m 1 Xffm 1 ,£ 2 m 2 (Q, 0)r° 2m2 y m i- (S24) 

l\m\i 2 ra 2 

Here we include only a static part v n = of xf^ti € 2 m 2 (^' ^«)' which corresponds to the clas- 
sical approximation. In this case, it is convenient to calculate not diagrams in Fig. S4a, but the 
maximum eigenvalue A of the Bethe-Salpheter equation (Fig. S4b), 

\4>tm{k) = V Uummi (k - p)T^G ei i>(p,iuj n )G m > mi (p + Qc,iuJ n )<Pe'm'(p), (S25) 

n 

where <pi m {k) = (flifk+Q c m) is the order parameter for the staggered pair, which has generally 
the momentum dependence. Indeed such /c-dependent order parameters have been proposed in 
some scenarios, including the unconventional spin density wave 8 , the orbital AFM 9 , the Helicity 
order 10 , and the spin nematic 11 . However, the present calculations indicate that 4>i m {k) with large 
eigenvalue A is almost k independent (not shown). Thus, it is appropriate that such order parameter 
is considered to be a kind of multipole, which is almost on-site pairing. The results for U = U' ~ 
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k+Q k'+Q ~- — < k+Q 

k+Q p+Q 




FIG. S4. Maki-Thompson type diagrams (a) and the corresponding Bethe-Salpheter equation (b). U — J 
phase diagram for JJ = U' + 2 J and J = J' at T = 16K (c). E~ and A 2 states fulfill A = 1 along red 
and blue lines, respectively. Broken lines are not realized. Darker background color corresponds to larger 
eigenvalue. Grey line denotes a change in the parameter set used in Fig. 3b. 

2.4 and J = J' = has been shown in Fig. 3a. Following Eq.(S22), /^-independent parts of the 
dominant two eigenstates, E~ and , are expanded as 

<j) E - = 0.150J x{y) - 0.21QT x(y)a + 0.206^(^)6 

+0.6711^! + 0.296D x{y)a2 + 0.592D x(y)b , (S26) 
4> A - = 0.825J, + 0.406T za - 0.140D 2al + 0.367D za2 . 

Thus the main ingredient of 4>e- is D x ^, while that of <\> A - is J z . This is the same as in the RPA 
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results. 

Fig. S4c is a U - J diagram for U — U' + 2J and J = J' at T = 16K. Large J stabilizes 
A 2 (J z ) state rather than E~(D X ^) state. Red and blue lines represent parameter lines fulfilled 
A = 1 for A 2 (J z ) and E~ (D x ( y )), respectively. Grey line corresponds to the parameter set used in 
Fig. 3b. 

Finally, within the same framework, we study the uniform magnetic susceptibility (Fig. 3c). 
This is obtained as Q = in Fig. S4a, 

IW> m >(£;) = n° m/ , m ,(A;) - U° £mAmi (k)V ilh;mim2 (k - k')U hm2;e > m >(k), (S27a) 
n°em,e>m>( k ) = -T^2G ee <(k,iu n )G mm <(k,iuj n ), (S27b) 

n 

X/.(0) = £ Jf^m,em'(k)J £ ; m '. (S27c) 
k 

B. Pseudospin representations 

When we focus on only ±5/2 orbitals, three multipole moments (D x , D y , J z ) defined in 
eqs.(S9), (S10) and (S20) are reduced to 

( 0.65 \ ( 0.65A ( 0.6 \ 

D x =\ \, D y = , J z =\ , (S28) 

^0.65 J y-0.65i J \ -0.6 y 

where these are normalized by eq.(S21). These are proportional to the Pauli matrices, {o x , a y , u z ). 
Thus we can describe these moments (D x , D y , J z ) as pseudospins approximately. 



VI. FERMI SURFACES IN ORDERED STATES 



Here we discuss the FS topology in some possible multipole ordered state. The FS was obtained 
by applying a finite effective field of the corresponding multipole. Figure S5a depicts the FS in the 
AFM order (A 2 ). We set 0.024 , I' 1 as the effective field. In this case, as shown in the left figure 
of Fig. S6a, the opening gap is of the order of 4meV, considering the renormalization of 1/10. 
The FS is similar to that in the previous DFT calculations 312 , except for the cage around T, which 
vanishes for larger effective field. Figure S5b represents the FS in the E~ dotriacontapole state, 
which is the same one as that in Fig. 4a. We set 0.024 (\1>|_ + \1>|_)/a/2 as the effective field, 
corresponding to the experimentally observed 'nematicity' along [110]. To clarify differences in 
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the FS topology, we use the same magnitude of effective field as the above A 2 AFM order. The 
FS in the E~ state is similar to that in the A 2 dipole state except for small separation of two FSs 
around M and the in-plane four-fold symmetry breaking of the cage FS, which is remarkable in 
the comparison of £ line (blue line) and E' line (green line) in the left figure of Fig. S6b. Figure 
S5c is the FS for the A{ state with 0.024\&^_. The FS around M is almost equivalent to the FS 
around X in the paramagnetic state. Finally, in Fig. S5d, we show the FS in the E + state with 
0.024(^| + + V% + )/y/2. In this case, the separation between two FSs around M becomes larger. 
Thus, the sizable change in different ordered states may be observed in detailed analysis of the FS 
topology. 
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FIG. S5. Fermi surfaces in ordered states, colored by the Fermi velocity. The right figures are two- 
dimensional cut for k z = 0. A2 dipole (a), E~~ dotriacontapole (b), A~[ dotriacontapole (c), and E + 
hexadecapole (d). The fourfold symmetry breaking in E~ and E + states is verified from anisotropy of the 
cage FSs in (b) and (d). All FSs have only slight differences. The most remarkable is a change of separation 
between two electron FSs around M. 
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FIG. S6. Dispersion relations along high-symmetry line in ordered states corresponding to Fig. S5. Red 
(blue) line represents dispersion relation in the paramagnetic (ordered) state. The energy scale is reduced 
by a factor of 10 as compared with Fig. SI a, by taking into account the mass renormalization effect. Left 
figures are the enlarged figures around the band crossing along T — M line. A large gap (~ 4meV) opens 
as compared with the paramagnetic band (red). Blue (green) line is a dispersion relation along E (£') line. 
In the E states, £ line and £' line are not equivalent. 
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